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Abstract. A characteristic feature of rainfall statistics is that they depend 
on the space and time scales over which rain data are averaged. A previously 
developed spectral model of rain statistics that is designed to capture this 
property, predicts power law scaling behavior for the second moment statistics 
of area-averaged rain rate on the averaging length scale L as L — > 0. In the 
present work a more efficient method of estimating the model parameters is 
presented, and used to fit the model to the statistics of area-averaged rain 
rate derived from gridded radar precipitation data from TOGA COARE. 
Statistical properties of the data and the model predictions are compared 
over a wide range of averaging scales. An extension of the spectral model 
scaling relations to describe the dependence of the average fraction of grid 
boxes within an area containing nonzero rain (the “rainy area fraction” ) on 
the grid scale L is also explored. 


1. Introduction 

Rainfall is a complex phenomenon involving the in- 
terplay of many physical processes in the atmosphere. 
Consequently, one often has to resort to a probabilistic 
description of the overall space-time distribution of pre- 
cipitation rather than a detailed physical model of the 
underlying processes. In contrast to physical models, 
statistical models have the advantage of being concep- 
tually economical in the sense that they generally in- 
volve fewer adjustable parameters. Moreover, they can 
be easily validated by comparing against data gathered 
over a large area and a large time period representative 
of a certain climatological regime. Once the model pa- 
rameters are determined from a sufficiently large data 
set, the model provides an efficient method of describing 
various statistical properties of rainfall over areas with 
similar rain climatologies [see, e.g., Bell et al , 2001, 
hereafter referred to as BKK]. 

Rainfall, despite its apparently irregular nature, is 
correlated both in space and in time. Spatio-temporal 


variability of rain directly influences other hydrologi- 
cal processes, such as soil moisture transport and sur- 
face runoff. Statistical behavior of area- and/or time- 
averaged rain rate at different length and time scales 
can be an important diagnostic for comparing actual 
rain behavior with predictions from global climate mod- 
els. It has important consequences for the estimation 
of the so-called sampling error that arises in satellite 
measurements because of their intermittent and often 
incomplete coverage of a given area of the globe [Mc- 
Connell and North , 1987; Bell et al , 1990; Bell and 
Kundu , 1996, hereafter BK96]. 

In this paper we present a statistical model of space- 
time variability of rainfall in the mesoscale regime de- 
rived from ground based radar observations. In this 
model, Fourier modes of the rain intensity R(x, t) are 
treated as a set of random variables whose evolution 
is governed by a linear stochastic equation driven by 
a white noise source. This equation immediately leads 
to an analytical form of the rainfall power spectrum 
depending on just a few adjustable parameters. The 
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Fourier transform of the spectrum yields the space-time 
covariance function of the point rain rate field, namely 

c(x,t,x',t') = (R'(x,t)R'(x',t')) (1) 

with 

R'(x, t ) = R{x, t ) - (R(x, t )) , (2) 

and determines how rain statistics depend on the aver- 
aging length and time scales. Here the angle brackets 
represent a statistical average over a hypothetical col- 
lection or ensemble of rain fields with similar rain cli- 
matology. The model was used previously in BK96 in a 
study of the sampling error in preparation for the Trop- 
ical Rainfall Measuring Mission (TRMM), where it was 
quite successful in describing the Global Atmospheric 
Research Program (GARP) Atlantic Tropical Experi- 
ment (GATE) Phase I data. Here we test it with a grid- 
ded precipitation data set constructed from the radar 
images obtained during the Tropical Ocean Global At- 
mosphere - Coupled Ocean Atmosphere Response Ex- 
periment (TOGA-COARE) [Webster and Lukas , 1992]. 
This data set covers a longer period of observation than 
the GATE data set and has a finer spatial resolution. 
Since the data set provides “instantaneous” spatial av- 
erages on a regular rectangular cartesian grid analogous 
to the GATE data, it is well suited to studying the 
statistical behavior of the rain field averaged over var- 
ious spatial scales. A parameter estimation method is 
described in this paper that utilizes the analytical be- 
havior of the model and is much more convenient and 
robust than the method used in BK96 which involved 
considerable amounts of trial and error. 

There exists a large class of physical rainfall models 
based on random processes that may be called hier- 
archical cluster models. These models attempt to cap- 
ture the spatial and temporal organization of storm sys- 
tems that cause precipitation. A general mathematical 
framework for such models was laid by Le Cam [1961] 
in which stochastic Poisson processes with suitably ran- 
domized parameters are used to describe the irregular 
occurrence of rain events in space and time. There is an 
extensive literature on this type of model [e.g., Gupta 
and Waymire , 1979; Cox and Isham , 1988; Smith and 
Krajewski , 1987; Waymire et a/., 1984]. However, these 
models generally involve a large number of parameters 
whose estimation from actual rain data presents a dif- 
ficult problem. For a recent account of a number of 
models in this category, including the problem of pa- 
rameter estimation in the context of a precipitation data 
set obtained in the United Kingdom from the HYREX 
project, a large network of radar and rain gauges, see 


Wheater et al [2000]. Because of the ease of param- 
eter estimation for real rain data sets, combined with 
its ability to capture the variability on many different 
space and time scales, we chose to consider a rainfall 
model based on a simple stochastic dynamical equation 
rather than a cluster model of the type described above. 

Other stochastic dynamical rainfall models exist in 
the literature [North and Nakamoto , 1989; Yoo et al , 
1996] which also describe the GATE spectrum with 
varying degrees of success. Of these the North-Nakamoto 
(NN) forced diffusion model is a special case of the 
model considered here. It describes precipitation as the 
diffusive spreading of rain events through a turbulent 
atmosphere. The Yoo- Valdes-North (YVN) model adds 
a second time derivative term which may help in model- 
ing the high frequency behavior of the rainfall spectrum. 
However, the model has more parameters and has so far 
been applied only in the low wave number regime, and 
it is not yet clear how well the model describes rainfall 
statistics at the smaller spatial scales considered in this 
paper. 

The present model is an outgrowth of previous work 
by Bell [1987] and Bell et al [1990] who constructed 
a model to simulate the statistics of the gridded rain- 
fall data from GATE. They introduced a gaussian ran- 
dom field on a discrete 4 km grid whose Fourier am- 
plitudes were assumed to obey a stochastic dynamical 
equation similar to the one used here. A nonlinear map- 
ping was used to generate the gridded rain rate field 
with a rain-rate probability distribution matching the 
observed lognormal distribution. Unlike that work, the 
present model directly describes the stochastic evolu- 
tion of the point rain rate field, which can be aver- 
aged to any desired length and time scale. Moreover, it 
yields the full space-time covariance function in terms 
of a small number of model parameters. Once they are 
determined from data, the model allows us to relate 
statistics at different scales, whereas the original model 
of Bell [1987] was applicable only to a single spatial grid 
scale and, moreover, the form of the spatial covariance 
function was specified as an explicit function of spatial 
separation, in effect increasing the number of empirical 
parameters. However the present model pays a price for 
its simplicity in being restricted to describing only the 
second moment statistics of rain. 

Our model spectrum naturally leads to power law 
scaling behavior of the rain statistics at small spatial 
scales suggesting an underlying fractal or multifractal 
structure of the rain field. Multifractal rainfall mod- 
els have been proposed by Schertzer and Lovejoy [1987] 
and Gupta and Waymire [1990] extending earlier ideas 
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of Lovejoy [1982], Lovejoy and Mandelbrot [1985] and 
Waymire [1985]. These models describe the moments of 
a spatially averaged rain field in terms of multifractals 
generated by an underlying self-similar multiplicative 
random cascade process [Marshak et al , 1994; Menabde 
et al, 1997]. The notion that the rain field is scale 
invariant over a broad range of spatial and temporal 
scales is inspired by statistical theory of turbulence and 
is supported by a growing body of empirical evidence 
[Foufoula-Georgiou and Krajewski , 1995]. An intrigu- 
ing feature of the multiplicative cascade models is that 
they naturally lead to a lognormally distributed precip- 
itation field. More recently, the spatial cascade models 
have been extended to space-time models which incor- 
porate dynamics in the scaling regime, either through 
a combined scaling of the space and time variables 
[Marsan et al , 1996] or through a description of the 
time evolution of the cascade generators in terms of 
an autoregressive process [Over and Gupta , 1996; Seed 
et al, 1999]. The multifractal approach has been ap- 
plied by Deidda [2000] to the GATE data to disag- 
gregate the large scale rain and compare the observed 
and synthetically downscaled rain statistics at smaller 
scales. Empirical evidence of scaling under a combined 
space-time scale transformation (“dynamical scaling”) 
has been presented by Venugopal et al [1999]. A review 
of the multifractal cascade dynamics models is given in 
[Schertzer et al, 1997]. 

The paper is organized as follows: Section 2 de- 
scribes the rainfall model. Here we collect a number of 
mathematical formulae related to the statistics of area- 
averaged rain rate fields that are useful for fitting the 
model to data. Of particular interest to us is the way 
these statistics depend on the averaging length scale. 
Readers who are mainly interested in the applications 
of the model to data may wish to omit the theoreti- 
cal details of the model and proceed to section 3 where 
we describe the gridded TOGA-COARE precipitation 
data set and the procedure for fitting the model to data 
and estimating the model parameters, returning later 
to section 2 for the main formulae as needed. Section 4 
is devoted to a discussion of the results, focusing in par- 
ticular on the scaling behavior of rain statistics at small 
spatial scales perhaps associated with fractal nature of 
the rain rate field. Some concluding remarks and direc- 
tions for future work are presented in the final section. 
Three appendices provide mathematical details of the 
derivation of a number of results utilized in the main 
text. 


2. Modeling Rainfall Variability 

2.1. A Spectral Model of Rain Rate Covariance 

In this subsection we give a brief account of the 
model. A set of results that are needed in fitting the 
model to data and checking some of its predictions are 
derived here. The underlying assumptions of the model 
as well as many of its salient features have already been 
described in detail in BK96. 

a. Model Basics 

As was emphasized in Bell [1987] and Bell et al. 
[1990], rain rate fields share with turbulent fluids the 
characteristic feature that the time scale over which spa- 
tial averages remain correlated depends on the size of 
the averaging area: the larger the averaging area, the 
longer the time over which the average rain rate re- 
mains correlated. This behavior must be captured by a 
realistic model of rainfall statistics. Although local rain 
statistics are expected to be affected by inhomogeneous 
and nonstationary directional effects such as occur, for 
instance, in rain bands, the model applications we en- 
visage require a description that is valid over a large 
enough region and sufficiently long time for such direc- 
tional effects to average out, so that the rain statistics 
can be treated to a good approximation as spatially ho- 
mogeneous and isotropic and stationary in time. These 
assumptions imply that the space-time covariance func- 
tion of the rain field at the points x, x' and times t,t f 
has the form 

c(x, t; x', t') = c(p, t) = c(p, t) 

where r = t' - t, p = x' - x and p — |p|. It can be ex- 
pressed as a Fourier integral 


C(P ’ T) = (2^372 I d 2 ke^P-^c(k,uj) , (3) 

where c(k,o;) represents the (unnormalized) rain rate 
power spectrum. Isotropy implies c(k, uj) = c(k,u) 
where k — |k|. The spatial Fourier amplitudes 

a(lc, t) = d 2 xe- ikx K(x,t) 

of the rain field fluctuation R'(x, t) defined by Eq. (2) 
are assumed to obey the first order linear stochastic 
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equation 1 


JU(M) = --a(k,t) + /(k,t). (4) 

at T fc 

[One could add a term (zk.v)a(k,f) to Eq. (4) to take 
into account the effect of advection, at the expense of 
isotropy.] Eq. (4) represents what is sometimes referred 
to as a first order autoregressive [AR(1)] process [Jenk- 
ins and Watts , 1968]. It is formally analogous to the 
familiar Langevin equation for the velocity of a particle 
executing thermal Brownian motion in a viscous fluid. 
Here /(k,£) is a white noise force satisfying 

</(k, t)f*( k', t')) = (2tt)F 0 8( k' - k) 5(t' - 1) , (5) 

where in the right hand side we have introduced the 
Dirac 6-function, and the asterisk denotes complex con- 
jugation. In Eq. (4) represents the characteristic re- 
laxation time of a(k, t). As in BK96, it is assumed to 
have the form 


Tfc " (1 + fc 2< L§) 1+1 ' ’ (6) 

where To and Lo are characteristic time and length 
scales. The length scale Lo effectively separates the 
long and the short wavelength regimes of the spatial 
fluctuations of the rain rate field. As a consequence of 
Eqs. (4) and (5) the power spectrum has the typical 
“red noise” form 


c(fc,u;) 


Fp 

^ 2 + ^“ 2 ‘ 


(7) 


The analytic form of r* in (6) is suggested by the form 
chosen by Bell [1987] and Bell et al. [1990] who, in 
turn, were motivated by the empirical results of Laugh - 
lin [1981] from GATE data. As already noted, it con- 
tains the North-Nakamoto model as a special case with 
v — 0. It may be pointed out that in the original work 
of Bell [1987] and Bell et al [1990], the form of the 
spatial covariance function was separately introduced 
from an empirical fit, whereas in the present model it 
is implicitly contained in the analytic form of r^. 

Because of the fractional exponent v, the implied 
stochastic dynamical equation governing the evolution 
of the rain rate field F(x, t) now takes the form of an 


integro-differential equation with an associated nonlo- 
cal character rather than a pure differential equation 
as in the case of NN and YVN models. It can also be 
formally expressed in terms of fractional order deriva- 
tives [see, e.g., Benson et al ., 2000]. Physically it can 
be thought of as representing a generalized diffusion 
process by which rain moves through a turbulent atmo- 
sphere. In particular, the negative values of v required 
to fit the data sets studied here appear to better de- 
scribe the high variability of rain rates at small spatial 
scales - an important factor in extrapolating statistics 
to rain gauge scales. 

The four model parameters Fo, To, L 0 and v char- 
acterize the second moment statistics. No specific as- 
sumption about the actual probability distribution for 
the random variable F(x,£) is introduced here. It 
should be noted however that the model of rain statis- 
tics remains incomplete without such an assumption, 
since the higher order moments cannot be represented 
within the model. Observationally, area averaged rain 
rates seems to be approximately lognormally or gamma 
distributed, at least on smaller space and time scales. 

b. Spatial Covariance Function 

Substituting the explicit form of the model spectrum 
in (3) and performing the integral over u> and the an- 
gular integral in the k-plane yields 

rOO 

c(p, t) == / dkkJ 0 (kp)c(k,T) , (8) 

Jo 

where Jo(x) is the Bessel function of order zero, and 

c(fc, t ) = y/^/2 F 0 Tfce _|r|/rfc (9) 

represents the lagged covariance of a(k, t) defined by 
the equation 

(a(k, £)a*(k',t')) = (27r)c(fc, t) <5(k' - k) , (10) 

and is the inverse temporal Fourier transform of c(fc, lo). 
At zero lag (t = 0), carrying out the remaining k- 
integration in (8), the spatial covariance of the point 
rain rate can be expressed in the simple closed form 

c(p, T = 0) = 70 C u {p/Lo) , (11) 

where, for convenience, 70 is defined by the relation 


1 A more general possibility is an equation of the form 

JU(k, t) = - [ d 2 k' — a(k', t) + / (k, t) , 

at J T kk > 

which allows coupling among different Fourier modes. 


F 0 = y/2fiT{l + v){LHt 0 ) 70 , (12) 

r(z) = f^°t z ~ 1 e~ t dt is the Euler T-function, and 

C v (z) - (z/2) t/ K u (z) (13) 
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is a function related to the modified Bessel function 
K„(z). For z«l, C v (z) has the behavior 

C v (z) = 5 r (^) + 2 T ^ (2) + ° {z2+2V) + ° (z2) ■ 

(14) 

For v > 0, the first term dominates and C v {z ) — > ^T(u) 
as z — ► 0, and the variance of the point rain rate 
c(0, 0) = ^7oI» is finite. However, as we shall see, 
like the GATE data, the TOGA-COAEE data requires 
v < 0, and the covariance function exhibits a power law 
singularity 

c(p,0)~i 7 or(|H)(p/2io)- 211 ' 1 (15) 

as the separation p — > 0. This remarkable singular be- 
havior can be traced to the fc-dependence of the char- 
acteristic time scale 77 governing the decay rate of the 
Fourier mode a(k, t) : For the long wavelength modes 
(k « 0) Tk approaches a constant value To, whereas for 
the short wavelength modes ( k » 1/To) r k behaves like 
£-2(1-11/1). Note that for values of the exponent v in 
the range — 1 < v < 0 (as seems to be the case for 
both GATE and TOGA-COARE data sets), -> 0 as 
k —> 00. For the NN model (u = 0) the singularity in 
c(p, 0) weakens to become logarithmic. 

Although the point variance is infinite in the model, 
this is not necessarily inconsistent with experience since 
observations always involve smoothing of the point 
rain rate field in space and/or time, and variances of 
smoothed rain field are finite in the model. 

c. Statistics of Area-Averaged Rain Fields 


Next, we turn to the space-time covariance function 
of rain rate averaged over a square area A = L 2 , 

7^4 (t) = -7 [ d 2 xi?(x,t) 

A J a 

at lag t namely, 


<Ua'(s,t) = (ft^(t)7^|,(t + r)) 

where A, A! are two L x L boxes whose centers axe 
separated by s and the prime on the rain rate variables 
denotes deviation from the mean, as in Eq. (2). Starting 
from the formula 

Caa'(s, t) = J^J^d 2 x J d 2 x' c(s + x' - x, r) (16) 
the spectral model yields at zero lag the expression 

C>tX'(s,0) = 70 J J d£id£ 2 (l - l£il)(l ~ 161 ) 


x C v 


* + Z 


F) . 

Lo) 


(17) 


See Appendix A in BK96 for details of the derivation. 
Upon setting s — 0 in Eq. (16), the variance of the area- 
averaged rain rate in a grid box of area A = L 2 is given 
by 

= <7 2 (L)=4 7 oG(^L/Lo), (18) 

where G(v; z ) is the double integral 

G(ir,z) = f 1 fdZidtii 1-6X1 -6) 

Jo Jo 


X C, 




(19) 


The power law singularity of the spatial covariance of 
the point rain rate as a function of separation leads to a 
simple scaling behavior of the statistics of area averaged 
rain with the size of the area. Using the first two terms 
of the expansion (14) in Eq. (19), which dominate as 
long as —1 < v < 0, we have 

G(^;z)«ir(-M) + ir(H)fr( 0 (|)" 2M , (20) 

where H(u) denotes the integral 

H{v)= f / 1 #id&(l-6)(l-k)(tf + & 2 )‘'- (21) 
Jo Jo 

Hence, for small L (more precisely L L 0 ), <r 2 (T) can 
be expressed in the form 

<t 2 (L) « a 0 + 6 0 T~ 2|1/I , (22) 


where ao and bo are constants given by 

ao = ^ 7o r(-M) , (23) 

60 = 2 7o r(M)tf (i/)( 2 To) 2M • (24) 

Thus the model predicts that as L — > 0, the variance 
of area-averaged rain rate <t 2 (L) exhibits a power law 
singularity: 

a 2 (L) ~ L~ 2 M . 


d. Correlation Length and Time Scales 


The lagged correlation function 


$AA'(s,r) 


Ca4'(s,t) 


a 


2 

A 


(25) 
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is a rapidly decreasing function of the separation s. The 
model correlation function is markedly nonexponential. 
Its fall-off rate is conveniently characterized by the “in- 
tegral correlation length” A int (L) [which would be the 
(1 /enfolding distance for an exponentially decaying cor- 
relation function] defined in the isotropic case through 
the relation 

A? nt (L) = rdss$ AA ,(s,T = 0) . (26) 

Jo 

It is related to the model spectrum through 

A? " t(L) = ^Z) c(k = 0 ’ r = 0) ’ (27) 

which upon using Eq. (18) reduces to 

A? nt (L) = ir(l + v)L\jG{y\ L/L 0 ) . ( 28 ) 


Prom Eq. (20) it then follows that Ai nt (L) vanishes as 
l ) v I as L — 0. 

The case v = — 1 presents a particularly simple and 
analytically tractable special case of our spectral model. 
In this case the spatial covariance of the point rain rate 
field satisfies c(p,0) oc 6(p). As a consequence, the cor- 
relation between the rain rates in disjoint grid boxes 
vanish. The integral correlation length Ai nt (L) is sim- 
ply equal to the box size L and the variance of area- 
averaged rain rate cr 2 {L) is proportional to 1/L 2 . 

We can also define an integral autocorrelation time 
Tint{L) (which would be the 1/e-folding time for an ex- 
ponentially decreasing autocorrelation function) through 
the formula 



dr <S> aa (s = 0,t) . 


(29) 


2.2. Extension to Other Rain Statistics 

The spectral model presented above makes a number 
of predictions about the change of rain statistics with 
scale. With a few additional assumptions these predic- 
tions can be shown to have interesting consequences for 
the asymptotic behavior of the “conditional” statistics 
of rain as well. We now turn to a brief description of 
these statistics. 

Consider remote sensing measurements of precipita- 
tion over a given area and time period and the resulting 
rain maps gridded on an L x L square grid. The ensem- 
ble mean rain rate ( R } is estimated as the average of all 
the data within the data set, and is independent of the 
resolution scale L . However, other statistical properties 
of the field, such as the fraction of grid boxes containing 
nonzero rain (“the rainy area fraction”) p(L) and the 
variance of the grid-box-averaged rain rate cr 2 (L), al- 
ready considered, depend nontrivially on L. The quan- 
tity p(L) representing the probability of nonzero rain 
in a grid box is a measure of the spatial intermittency 
of rain at the scale L. Two other related statistics are 
the mean r c (L) and the variance cr 2 (L) of the grid-box- 
averaged rain rate conditional on the box containing 
nonzero rain. These are related (see, e.g., BKK, Ap- 
pendix A) through 

(R)=p(L)r c (L), (32) 

a*(L) = p(L)[a 2 c (L) + r C (L) 2 } - ( R ) 2 . (33) 

It is convenient to define the dimensionless ratio 

p = p(L) = a c (L)/r c (L). (34) 

Using Eq. (32), Eq. (33) can then be written as 


It can be expressed in terms of the Fourier transform 
C A (u) of the lagged autocovariance function C AA ( s = 0, r) 
as 

WD = ■ (30) 

Evaluation of C A ( 0) (see Appendix A) yields 


Tint(L) = T 0 


T(1 4- v) G(1 + 2u\L/Lo) 
r(2 + 2z/) W^jLoT 


(31) 


The values of v that best describe rain data in fact ap- 
pear to lie in the narrower range — 1/2 < i/ < 0 (rather 
than the range — 1 < v < 0 considered above). The 
numerator of Eq. (31) is then finite and consequently 
the model predicts that r int (i) — > 0 as L 2 ^ in the Unit 
L — ► 0. In other words, the time scale for fluctuations 
of average rain rate in a box of area L 2 gets smaller 
with decreasing box size. 


<t 2 (L) 1 +£ 

(R) 2 P(L) ‘ 


Based on analysis of rain gauge data from Darwin, 
Australia and Melbourne, Florida and radar data from 
GATE, Short et al [1993] have suggested that p is rela- 
tively constant over a range of length scales, averaging 
times, types of data (radar or rain gauge) and rain cli- 
mates. In particular, p is insensitive to the averaging 
length scale L. This was also found by BKK for SSM/I 
rain data. 

If the quantity p(L) remains finite as L — » 0, then 
Eq. (35) together with the asymptotic behavior (22) 
predicts that p(L) would vanish as L 2 ^ in this limit. 
The limiting behavior p(L) — ► 0 as L — ► 0 in turn im- 
plies that the instantaneous point rain rate field i?(x, t) 
vanishes almost everywhere , i.e., everywhere except a 
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set of measure zero - a familiar property of spatial frac- 
tals. We shall return to this in greater detail in section 
4 when we present the results of our data analysis. 

3. Testing the Model Predictions with 
Ground based Radar Data 

3.1. The TOGA-COARE dataset 

The radar data we use to test the spectral model 
described above were collected during TOGA- COARE 
[Webster and Lukas , 1992] carried out in the Western 
Pacific equatorial warm pool during the period Novem- 
ber 1992 to February 1993. The measurements were 
made with a pair of calibrated Doppler radars (labeled 
TOGA and MIT) on board two research ships located 
near the center of the Intensive Flux Array (IFA) - a re- 
gion in the shape of a quadrilateral, roughly centered at 
2° S, 156° E. Reliable rain rate estimates were provided 
out to a distance 145 km from each radar. The fields 
of view (FOVs) of the two radars overlapped substan- 
tially when they were both operating. The measured 
radar reflectances Z were converted into average sur- 
face rain rates on a 2 km spatial grid using an empirical 
Z — R relation. The data set is described by Short 
et al [1997]. The total 101-day period of observation 
was divided into three “monthly” cruises covering the 
periods November 11 to December 10 1992, December 
15 1992 to January 18 1993 and January 23 to Febru- 
ary 23 1993 designated respectively as Cruises 1, 2 and 
3. Rain maps were available roughly at intervals of ten 
minutes with occasional long gaps during which data 
from one or both of the radars were missing. 

3.2. Statistical Analysis of the Data 

The gridded data set obtained from the TOGA- 
COARE observations yields estimates of the various 
area-averaged rain statistics introduced previously, al- 
lowing one to test the model predictions. We now briefly 
describe how these statistics are obtained from the data. 

The statistics are derived using rain estimates from 
square areas 128 x 128 km 2 centered at the location 
of each of the two ships. The locations of the areas 
changed slowly with time because of the slight drift of 
the ships during the month-long experiments. However, 
the effect of this on the statistics, if any, is expected to 
be negligible. Each area is subdivided into boxes Lx L 
in size, with L — 2 n km (n = 1, 2, . . . , 7), and the 
quantities (R), a 2 (L ), p(L ), r c (L) and a c (L) defined 
in the previous section are computed for each L from 
the data. Only those boxes in a rain map were used for 


which at least 95% of the box had valid data. This was 
necessary since the radar images frequently had data 
missing along a radial line emanating from the center of 
the circular FOV where the radar was located. These 
dropouts affected the grid boxes located close to the 
origin through which the radial line passed, especially 
at the smaller grid length scales L = 2, 4, 8 and 16 km. 

For comparison purposes, we also computed the same 
spatial statistics for the GATE Phase I (1716 images) 
and Phase II (1512 images) data from the 280 km box 
centered at the location of the ship. For GATE data 
the basic grid size is 4 km. 

We determined the integral correlation time T- mt (L) 
for each L as follows. Mean rain rate time series were 
constructed for each Lx L box within the selected areas 
for the three cruises. Again only boxes with at least 95% 
of the area containing data were included in the time se- 
ries in order to contend with the data gaps mentioned 
above. The lagged autocovariance function Caa(0, t) 
averaged over all boxes in the chosen areas was com- 
puted for each L at lags r of multiples of 10 min. 

As an example, the lagged autocorrelation function 

* L ( r ) = $aa(0,t) 

for L = 128 km for the six data sets are shown in 
Figure 1. For each data set a numerical estimate of the 
time integral T int (T) = / 0 °° $L(r)dT was obtained. This 
is discussed in more detail later. 

The next subsection describes how the parameters 
7o, Lq, v and tq are determined by fitting the data to the 
statistics of gridded averages predicted by the model. 

3.3. Estimation of the Model Parameters 

The parameters 70 , To and v determine the spatial 
statistics at zero lag. They are estimated by fitting the 
variance cr 2 (L) to the asymptotic formula (22) valid for 
small L. The fit determines the exponent v and the 
coefficients ao and bo. Eq. (23) determines the parame- 
ter 70 and Eq. (24) then determines Lq. The procedure 
was carried out for the TOGA-COARE as well as both 
GATE Phase I and Phase II data sets. 

Once v and Lq are determined, an estimate of the 
parameter tq for each data set can be obtained by aver- 
aging over the values individually computed by solving 
(31) for To using the values of Ti nt (T) for various box 
sizes L. However, since there is only a finite length T 
of data available, instead of ri nt (T) one can only obtain 
an estimate 

fTm ax 

Tmt(I',T 0 ;T max ) = / $x,(r)dT, 

Jo 


( 36 ) 
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Lagged Autocorrelation for TOGA CO ARE Data Sets 



Figure 1. Plot of lagged autocorrelation functions 
$l(t) for L = 128 km square boxes centered at the 
location of the TOGA and MIT radars for each of the 
three cruises. 

with a finite cut-off r max < T imposed on (29). In prin- 
ciple, we could solve (36) for r 0 using the model from 
the data-derived value of fi nt for a chosen cut-off r max . 
In practice this is difficult, since the model-predicted 
7o-dependence of the function fi nt (T,To;T max ) is non- 
linear and complex to calculate. Instead, we adopt an 
iterative method described next. 

There is an uncertainty inherent in choosing r max in 
order to get the “best estimate” of To. It arises because 
the tail of the autocorrelation function suffers from a 
combined effect of statistical noise and possible secular 
variations of the rainfall pattern over longer time scales 
not represented in our spectral model. The choice of the 
cut-off is guided by several considerations. On the one 
hand, r max should be large enough so that the contribu- 
tion of the tail in the integral (29) is small. Some guid- 
ance on an acceptable value of T max is provided by the 
model’s prediction of this contribution. However, the 
empirical lagged autocorrelation suffers from sampling 
error because it is obtained from a finite time series of 
length T « 30 days. The cut-off T max must therefore 
satisfy the condition T max T. Also r max should be 
chosen to be small enough so that the effect of slow sec- 
ular modulations of the lagged autocorrelation function 
caused by physical phenomena such as a diurnal cycle 
or various atmospheric waves, which are not accounted 
for in our simple model spectrum, is minimized. Evi- 
dence of such longer time scales was most obvious in the 


Cruise 2 MIT and the Cruise 3 TOGA radar data sets. 
These modulations can probably be attributed to var- 
ious large scale atmospheric waves known to influence 
rainfall in the equatorial Pacific. We shall revisit some 
of these issues in the next section. Keeping these factors 
in mind, we chose a relatively small cut-off T max = 12 
hours with the implicit assumption that the spectral 
model accurately describes the observed correlations up 
to the lag r = r max . 

Next, we turn to the problem of estimation of To. 
Using (36), Eq. (29) for the integral correlation time 
^int(T) = Ti n t(L,To) can be expressed as the sum 

poo 

7int(T, To) = Tint(T,To; Tmax ) + / $l(,t)cIt (37) 

** T'max 

If the left-hand side of (37) were known, Eq. (31) would 
immediately yield To. However, while the first term 
on the right-hand side of (37) is obtained from data, 
the second term can be calculated from the model only 
when To is known. We therefore solve Eq. (37) itera- 
tively: to begin the process, we obtain a first guess for 
To by ignoring the contribution of the integral; for all 
the following steps of the iteration we evaluate the inte- 
gral on the right-hand side using the model $j,(t) with 
the value of t 0 obtained from the previous step. For 
ease of computation, we replace the square box by a 
circular disk of equal area. This approximation proved 
to be accurate enough for our purposes. Appendix B 
gives some details of the model calculations for a circu- 
lar area. In practice, for the TOGA-COARE data, the 
correction to r 0 due to the introduction of a cutoff is of 
order 5%, and it was unnecessary to go beyond the first 
iteration. 

The iterative method applied to GATE Phase I data 
with r max chosen to be 18 h gives an estimate of t 0 
within about 1% of what was reported in BK96 from a 
direct fit to the lagged autocorrelation function. 

4. Results and Discussion 

The full set of parameter values for which the model 
best fits the TOGA-COARE data are given in Table 1, 
along with the corresponding values for GATE Phase I 
determined in BK96 using a somewhat less exact and 
more arduous trial-and-error curve fitting method in- 
stead of the analytical method adopted here. 

4.1- Model Parameters v , L 0 and 70 

Table 1 summarizes the TOGA-COARE model pa- 
rameters. We see that the values of the power law expo- 
nent v range between about —0.21 to —0.34 indicating 
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Table 1. Parameter Values for Rain-Rate Covariance Model. 


Dataset 

7o (mm 2 h -2 ) 

V 

Lq (km) 

T 0 (h) 

(R) (mmh x ) 

GATE Phase I 

1.0 

-0.11 

104. 

13.0 

0.492 

TOGA Cruise 1 

0.067 

-0.335 

94.06 

6.8 

0.139 

MIT Cruise 1 

0.086 

-0.297 

73.89 

5.8 

0.134 

TOGA Cruise 2 

0.616 

-0.239 

53.81 

5.0 

0.352 

MIT Cruise 2 

0.206 

-0.205 

70.40 

8.2 

0.229 

TOGA Cruise 3 

0.127 

-0.290 

61.04 

5.2 

0.155 

MIT Cruise 3 

0.180 

-0.259 

64.94 

4.5 

0.200 


Parameters for the model spectrum obtained from fits to radar data from Phase I of GATE and from the 
ships MIT and TOGA during the three TOGA-COARE cruises. Average rain rate for each dataset is given in 
the last column. 


a much stronger singular behavior than was found for 
the GATE Phase I data (v = —0.11). The characteristic 
time To lies in the range of roughly 4.5 to 8.2 hours, con- 
siderably shorter than the GATE Phase I value To « 13 
hours. The characteristic length Lq ranges between 
about 54 and 94 km indicating a somewhat sharper 
fall-off of the spatial correlations compared to GATE 
Phase I for which Lq = 104 km. The results exhibit- 
ing the quality of the fit are shown in Figure 2 for the 
TOGA and MIT radars. It is seen that in most cases 
the fit is quite good up to grid size L = 64 km, beyond 
which the approximation on which Eq. (22) is based is 
expected to break down. Our results seem to be consis- 
tent with the model prediction that the variance <r 2 (L) 
has a power law behavior a 2 (L) ~ L~ 2 ^ as L -* 0. Ro- 
bustness of the parameterization was checked by com- 
paring the variance of area-averaged rain rate (o^th 
predicted by the exact formula (18) for an area A = L 2 
with L = 128 km against the observed variance (cr^) 0 bs 
for the 128 km box centered at the two ship locations. 
Table 2 shows the comparisons, which are quite good. 

For GATE data we find that for L between 4 km and 
56 km, the variance cr 2 (L) is quite accurately fitted (not 
shown here) by the formulas 

2 ( t\ f 17.3L -0,24 — 4.66 mm 2 h' 2 , Phase I, 

a K ’ \ 15.0 £-°- 38 - 1.92 mm 2 h- 2 , Phase II. 

Equations (22) - (24) yield the parameter values 

7o = 1.02 mm 2 h" 2 , v — —0.12, Lq = 93.8 km 

for Phase I and 

7o = 0.63 mm 2 h" 2 , v = -0.19, Lq = 82.7 km 


for Phase II. [The original GATE Phase I parameters 
listed in Table 1 lead to the asymptotic formula (t 2 {L) = 
16.80L' 0 * 22 — 4.89 mm 2 h" 2 (see BK96 Appendix A), 
whose coefficients and the power law exponent are well 
within the 95% confidence interval of the corresponding 
quantities in the empirical least squares fit given above.] 

We also explored the dependence of the model pa- 
rameters on the mean rain rate (R). While no system- 
atic dependence could be discerned for Lq , to and 17 
the strength parameter 70 was found to depend on ( R ) 
according to the simple formula 

70 = <R)\ (38) 

where k is a dimensionless constant. The fit in Figure 
3, which also includes points corresponding to GATE 
Phase I «i?) = 0.492 mmh' 1 ) and GATE Phase II 
(( R ) — 0.386 mmh -1 ), yields a value k « 4.33. 

4.2. Model Parameter t 0 

Next, we turn to a discussion of the dependence of 
the integral correlation time T[ n t(L) on the spatial av- 
eraging scale L. The plots in Fig. 4 show reason- 
able agreement between the values of Tj nt (L) determined 
from data and the corresponding theoretical predictions 
computed for various box sizes L using the model pa- 
rameters. Linearity of the log-log plot for small L con- 
firms the power law dependence T mt (L) ~ ( L/Lq ) 2 ^ 
as L — > 0 expected from our model. In BK96 it was 
found that the integral correlation time estimates for 
the GATE Phase I is described quite well by the power 
law form r [nt (L) = 0.67 L 0A9 h. Again the exponent of 
the power law is rather different from the expected value 
0.22 obtained from the spectral model fit. However, it 



10 


Table 2 . Comparison between theoretical and observed values of g\ for L = 128 km. 


Dataset 

(oi)tfa (mm 2 h 2 ) 

(^)obs (mm 2 h 2 ) 

TOGA Cruise 1 

0.107 

0.116 

MIT Cruise 1 

0.093 

0.105 

TOGA Cruise 2 

0.399 

0.371 

MIT Cruise 2 

0.176 

0.174 

TOGA Cruise 3 

0.107 

0.119 

MIT Cruise 3 

0.155 

0.169 


Comparison between values of the variance of area-averaged rain rate g\ computed from the model using 
Eq. (18) and the parameters listed in Table 1 for L = 128 km and from the data for two 128 km boxes centered 
at each of the two ship locations. 


was also found there that while the model captures the 
overall decay of the autocorrelation function it 

substantially underestimates the correlation at smaller 
lags, indicating a departure from the model spectrum 
for a first order process. Indeed, it was shown by Bell 
and Reid [1993] that this behavior could be accounted 
for by going over to a second order model, which of 
course introduces more adjustable parameters. 

A number of caveats in the estimates of To should be 
pointed out. A fundamental uncertainty in the value of 
r^t (L) arises from the fact that it must be estimated 
from a time series of finite length T which in our case 
is of order 30 days. It is known [Bell, 1980] that if the 
true autocorrelation function of a normally distributed 
random variable has an exponentially decaying form 
$l(t) = exp(— r/rint) then the integral correlation fi nt , 
defined by Eq. (35) with a finite cut-off T max ? has vari- 
ance 

<7 2 [fint] = 4(r max /r)fj? t , (39) 

where it is assumed that T[ nt (L) r max T. A proof 
of Eq. (39) is given in Appendix C. Thus, in order to 
reduce the statistical uncertainty in the estimate, r max 
should be chosen to be as small as possible but large 
enough to include the entire range over which 4 >l(t) 
differs appreciably from zero. In our case, if we choose 
Tmax ~ 12 hrs, then it follows that the 2 <r uncertainty in 
the estimate of r\ ut (L) obtained from a single time series 
could, in fact, be as large as 0.5fi nt . It is to be noted 
that, strictly speaking, this estimate holds only for a 
gaussian random variable undergoing a linear stochas- 
tic dynamical process with an exponentially decaying 
autocorrelation function. However, the probability dis- 
tribution of rain rate variable is markedly nongaussian 


and the autocorrelation possesses a much more slowly 
decaying tail. These factors can lead to an uncertainty 
larger than what is predicted by Eq. (39). On the other 
hand, for the smaller grid boxes of size L = 2 n km with 
n < 6 , the actual uncertainty should be considerably 
lower than the estimate (39) since the autocovariances 
were obtained by averaging over the results for the time 
series for the 4 7-n boxes. These time series are, how- 
ever, not completely independent of one another, be- 
cause of the spatial correlation of the rain rate in a box 
with others in its vicinity, and an exact estimate of the 
uncertainty would be difficult to obtain. 

4.3. Conditional Statistics and the Spectral 
Model 

The spectral model described in section 2 allows one 
to quantify how the statistics of area-averaged rain rate 
depend on the averaging scale L . These statistical quan- 
tities, such as the variance <r 2 (L) and the integral cor- 
relation length and time scales, pertain to variability of 
the rain field including its zeroes, i.e., the region where 
it is not actually raining. However, unlike most random 
variables in nature, the instantaneous rain rate has a fi- 
nite nonzero probability of having the value zero, in ac- 
cordance with common experience - namely, that it is 
not raining in most places most of the time. Mathemati- 
cally this is expressed by the property that the probabil- 
ity density function for the area-averaged instantaneous 
rain rate TZa(^) defined by ( 1 ) has a ^-function at zero 
followed by an (approximately lognormal) continuum 
[Kedem et al , 1990], The probability p(L) for the value 
7 Za 7 ^ 0 is expected to have nontrivial L-dependence 
and is an important statistic of the gridded rain field. 
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Variance as function of area for TOGA COARE 




Figure 2. Comparison between the observed variance 
of area- averaged rain rate cr 2 (L) for the TOGA-COARE 
data and the corresponding spectral model predictions 
for various box sizes L ranging from 2 km to 128 km 
using Eqs. (22)-(24) and the parameter values listed in 
Table 1. 


1.2 

1.0 

0.8 


0.4 
0.2 
0.0 

0.1 0.2 0.3 0.4 0.5 0.6 

Mean Rain Rate <R> (mm h' 1 ) 

Figure 3. Plot of the parameter 70 against the mean 
rain rate (R) for the TOGA-COARE and GATE data 
sets. The continuous curve shows the least squares fit 
to the formula (38). 
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As discussed in section 2, an additional property of the 
conditional statistics of rain, namely the constancy of 
the dimensionless ratio /x, allows one to draw some infer- 
ences about the asymptotic L-dependence of the func- 
tion p(L). We now turn to a discussion of this statistic 
for the TOGA-COARE data set. It is, however, worth 
pointing out that in doing so we are exploring terri- 
tory that lies outside the scope of the original spectral 
model. 

Our data set provides a simple way to test the con- 
stancy of the ratio p defined by (34) as originally sug- 
gested by Short et al [1993]. Eq. (35) implies that the 
plot of 1 -1- o 2 (L)/(R) 2 vs. 1 /p(L) would be a straight 
line through the origin if p is a constant independent of 
L. A scatterplot of all the data points from the six data 
sets (2 ships, 3 cruises) and each value of L is shown 
in Figure 5. The linear fit seems satisfactory; the slope 
yields a mean value of p « 3.03 collectively for all the 
six data sets. However, there is a pronounced system- 
atic deviation from the linear fit in the “large” p(L) 
regime [i.e., p(L) ~ 1], which corresponds to large L. 
One can also evaluate p separately for each data set. 
The individual values were found to be within a narrow 
range, with a low value of 2.63 for MIT Cruise 2 and 
a high value of 3.39 for TOGA Cruise 2. As noted in 




Integral correlation time (h) Integral correlation time (h) 
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Plot of integral correlation time vs. L for TOGA COARE 




Box size L (km) 


Test of constancy of \i 



Figure 5. Scatterplot of 1 + a 2 (L)/(R) 2 vs. 1 /p(L) 
for all the six TOGA-COARE data sets for averaging 
length scales L ranging from 2 km to 128 km (6 x 7 = 42 
data points). Least squares fit to a straight line through 
the origin yields the slope 1 + fi 2 ~ 10.2, i.e. (i « 3.03 
for the entire TOGA-COARE data set. 


Figure 4. Plot of the integral correlation time Tint 
evaluated with a finite cut-off T max for the three TOGA 
and MIT cruises and various box sizes L ranging from 
2 km to 128 km. The curves represent interpolation of 
the values computed from the spectral model. 
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section 2, constancy of p (or more generally, finiteness 
of p as L — * 0) together with (22) leads to another im- 
portant conclusion - a power law scaling behavior of 
the rain probability p(L) at small L . Combining the 
formulas (22)-(24) and (38) with Eq. (35) we infer that 
for L <C Lo , 

l/p(i) oc 1 + /3(L/2L 0 ) _2|l/| ) (40) 


where 


2KT(\ V \)H(y) 

P l + KHH)’ 


(41) 


H(v) is given by Eq. (21) and k is the dimensionless 
number introduced in (38). This immediately implies 
the power law behavior 


p(L) ~ 



(42) 


as L — ► 0, which can be independently checked from the 
data. Fits to a power law p(L) oc L n for each of the six 
data sets are shown in Figure 6. A straight line fit on 
a plot of log p(L) vs. log L provides evidence of a power 
law behavior between L — 2 km and 64 km. In Table 
3 the observed exponent t] 0 ^ s is compared against the 
theoretically expected value 77 th = 2\v\ for each data 
set. We see that, although there is generally speaking a 
close agreement between the two sets of values, for four 
of the six data sets 77 ob s is somewhat larger than 77 th 
contrary to what one would infer from Eq. (35). How- 
ever, it should be noted that according to the model 
the power law dependence on L becomes accurate only 
in the scaling regime L <C Lo. Consequently, it is more 
appropriate to compare 77 th with the power law expo- 
nent determined from p(L) in the small L limit. Indeed, 
a closer inspection of the graphs of p(L) shows a slight 
but systematic deviation at small L from the power law 
obtained by least squares fit, indicating a smaller slope 
there. In Table 3 we have also included an estimated 
exponent r}® bs evaluated from the slope of a linear in- 
terpolation between L = 2 km and 8 km, which clearly 
improves agreement between the model and observa- 
tion. One can also see in retrospect that the constancy 
of p which gave rise to the power law scaling behavior 
of p(T), itself becomes progressively more accurate in 
the limit of small L. 

For the GATE data we find that within the range of 
L between 4 km and 56 km the rain probability p(L) 
(L in km) is represented quite well by the power law 
formulas 


Plot of p(L) vs. L for TOGA COARE 




Box size L (km) 


Figure 6. Power law fit to the fractional rainy area 
p(L ) for the three TOGA and MIT cruises and various 
box sizes L ranging from 2 km to 64 km. The exponents 
are given in Table 3. 


m = J 0.122(I/4)°- 61 Phase I 
P{ ) \ 0.098(L/4)°- 68 Phase II ’ 
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Table 3. Comparison between theoretical and observed values of the exponent rj . 


Dataset 

'Hobs 

^obs 

Vth = 2\u\ 

TOGA Cruise 1 

0.760 ± 0.026 

0.571 

0.669 ± 0.025 

MIT Cruise 1 

0.728 ± 0.025 

0.618 

0.595 ± 0.047 

TOGA Cruise 2 

0.464 ± 0.011 

0.382 

0.477 ± 0.057 

MIT Cruise 2 

0.540 ± 0.014 

0.417 

0.409 ± 0.054 

TOGA Cruise 3 

0.582 ± 0.015 

0.448 

0.581 ± 0.056 

MIT Cruise 3 

0.573 ± 0.015 

0.448 

0.519 ± 0.063 


Comparison between value of the power law index 7] defined by the stipulated behavior p(L) oc L 71 of the 
fractional rainy area p(L) as function of the averaging length scale L determined from the data and the value 
predicted by the model. The second column gives the “mean” exponent obtained by an overall least squares fit, 
while the third column gives the exponent determined from the slope at small L . In the last column specifying 
the “theoretical” value 7? t h = 2|i/| determined from the a 2 (L) vs. L plots [Figure 2 via the fit to (22)], we also 
include the error bars not reported in Table 1. 


which are consistent with the expected behavior p(L) — ► 

0 as L — ► 0. [See also Ha et al , 2002, who use a linear 
regression fit to describe p(L) as function of L]. The 
exponents are, however, in each case quite different from 
the values suggested by the spectral model fit. We also 
find that for GATE, the quantity (i{L) actually depends 
significantly on L . In fact, p(L) changes from about 
1.77 at L = 4 km to about 2.16 at L = 56 km for Phase 

1 and from about 1.89 at L = 4 km to about 2.22 at 
L = 56 km for Phase II. 

As in the TOGA-COARE case, the goodness of a 
simple power law fit is probably somewhat fortuitous 
for GATE as well, since deviations from a power law be- 
havior are to be expected on the basis of the model be- 
cause of the presence of “subdominant” terms in l/p(L) 
such as the constant term in (40) or other possible In- 
dependent terms arising from a weak L-dependence of 
p(L). Upon attempting to fit the data for p(L) to a 
form like (40) we find that the best fit values of the 
exponent r) are slightly smaller than the values given 
above. However, they still differ substantially from the 
model-inspired value 2\v\. We have been unable to fully 
reconcile the apparent discrepancy between the two sets 
of exponents along the lines discussed above in connec- 
tion with the TOGA-COARE data. 

One particular aspect that we did not address in this 
paper is the dependence of p(L) and other rain statistics 
on the value of the rain threshold, i.e., the smallest value 
of rain rate that is experimentally measurable. This is 
an important issue, since different instruments have dif- 
ferent thresholds and hence threshold dependence can 


greatly complicate intercomparison efforts for rain rate 
averages, especially conditional averages. 

4.4. Discussion 

From Table 1 it is evident that the model param- 
eters vary considerably among the data sets. These 
variations are almost certainly due in part to real sec- 
ular variations in the rainfall pattern. As discussed 
by Short et al. [1997], the TOGA-COARE region of 
the tropical Western Pacific experiences planetary scale 
equatorial waves including an annual cycle, intrasea- 
sonal oscillations such as Madden- Julian oscillations 
[Lau and Chan , 1985; Madden , 1986; Nakazawa , 1995] 
and westerly wind bursts related to El Nino-Southern 
Oscillation (ENSO) events. The first half of Cruise 2 
had an intensely active period associated with intrasea- 
sonal oscillation accompanied by passage of cloud su- 
perclusters and heavy rainfall. The rainfall time se- 
ries during this period also shows a quasi 2-day oscil- 
lation of the type described by Takayabu [1994] and 
Takayabu et al. [1996]. They attribute this oscillation 
to westward-propagating inertio-gravity waves. By con- 
trast, the spatial rainfall pattern during much of Cruise 
1 appears to be much more homogeneous and isotropic, 
characteristic of small-scale convective systems. For a 
more detailed discussion see Short et al. [1997]. 

5. Conclusion 

In this paper we have presented a model of second 
moment statistics of rainfall based on a linear stochas- 
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tic dynamical equation. The model makes definite pre- 
dictions about the statistical behavior of area-averaged 
rain rate in both space and time. In particular, the 
model predicts simple power law dependence of the rain 
statistics on the averaging area as the area tends to 
zero, perhaps reflecting a fractal nature of the rain field 
at smaller spatial scales. These predictions are tested 
with a radar-derived gridded tropical oceanic precipita- 
tion data set and the model is shown to describe the 
data quite well within the intrinsic limitations of the 
model assumption of statistical homogeneity, isotropy 
and stationarity. In particular, the intimate relation- 
ship between the dependence of the correlation length 
and time scales of the rain field on the spatial smooth- 
ing predicted by the model is in accordance with data. 
It should, however, be emphasized that in extrapolating 
the model-derived limiting behavior of the rain statis- 
tics deduced from a data set on a finite grid as the 
averaging scale goes to zero, we are effectively probing 
the “subgrid” structure of the rain field, i.e., the struc- 
ture on scales finer than allowed by the instrumental 
resolution. 

To conclude, we briefly mention some directions for 
future work. Like the spatially averaged rain rate, the 
spectral model also predicts relationships among the 
power law scaling exponents for the time-averaged rain 
rate as well. It would be interesting to see how well the 
model is capable of describing statistics of rain gauge 
data averaged over various time intervals. In a sequel to 
this paper, we will investigate the statistical properties 
of time-averaged local rain rates as measured by rain 
gauges using our model spectrum. Our model provides 
a convenient framework for addressing intercomparison 
problems that involve averaging the precipitation field 
over many different length and time scales. Recently the 
model has been used to investigate the problem of val- 
idation of satellite measurements of rainfall using rain 
gauge measurements on the ground [Bell and Kundu , 
2003] . Finally, simple power law behavior of the various 
statistics of rain fields with respect to rescaling of the 
spatial averaging grid at small scales suggests a model- 
independent explanation based entirely on the scaling 
properties of the various statistical quantities involved. 
We plan to explore this possibility in future publica- 
tions. 


Appendix A: Analytic Formula for 

Tint (^) 


In this appendix we give a derivation of the formula 
(31) for the integral correlation time Ti nt (L) and de- 


scribe its limiting behavior as L — ► 0. 

Equation (18) for the variance of the area- averaged 
rain rate in a grid box A of area A = L 2 reads 

a 2 A = 4^G(v-L/L 0 ), 


where G(v\z) is defined as the double integral (19). 
On the other hand, the lagged autocovariance function 
Caa{0,t) is given by (BK96, Appendix A) 

caa(o,t) = | J~J ~ dhdk2g2 (^r) 

xQ 2 (^r) c (k> r) 

where Q 2 (x) = sin 2 x/x 2 is the Bartlett filter function 
and c(k, r) is give by Eq. (11). Setting r = 0, we have 
an alternative integral representation of cr\ : 




' k 2 L\ 


FqTo j°°j°°dhdk 2 9 2 (hL/2)9 2 (k 2 L/2) 


[1 + (k 2 + kl)Ll) 


(Al) 


where in the last step we have made use of the explicit 
form of Tk given by (6). Comparing the two represen- 
tations of cr\ we obtain the useful integral identity 


n °° dk x dk 2 G 2 (fciL/2) Q 2 (hL/2) 
[i+w+tj>i§i>+- 

= mfu) TftoW- < A2 > 


Next, consider the Fourier transform of Caa( 0, t) : 

c~ah = i 

x9 2 c(k,u) . 

Using the explicit form of the model spectrum (7) at 
uj = 0, we get 

C a (lj - 0 ) = ^F 0 J^J^dkrdkz 9 2 (^) 
*9 2 (^) r 2 k 

2 2 f°° f°°dk 1 dk 2 g 2 (k 1 L/ 2 )g 2 (k 2 L/ 2 ) 

7T OT °Jo J 0 [l + (k 2 + k 2 2 )Lff(i+'') • 
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In view of the identity (A2) this reduces to 

C A (v = 0) = 4^-- 7o To p ^2 _^~2 ^ G ( 1 + 2 ^’ ’ 

(A3) 

where we have also used (12) to substitute for F 0 in 
terms of the standard model parameters 70, L 0 and v. 
Using Eqs. (18), (A3) and (30) we obtain Eq. (31). 
The integrals were numerically evaluated with the help 
of Mathematica [v. 4; see Wolfram , 1999]. 

Finally, we examine the limiting behavior of T in t(L) 
as L — > 0. We note that for the range of values of 
the exponent v of interest to us, — 1/2 < v < 0, (or 
0 < l+2i/ < 1), Eq. (20) yields the asymptotic behavior 

G{r,z) = ±T{\v\)H(v) (|)' 2M +0(1) 


The areal integrals can be done in a closed form. 
After some algebra using (9), one obtains 


C AA (0,r) 


4 7o r(i + v) r°°dK jI(kol) 

O 1 Jo K V ( K ) 

X e -\r/ro\v(K) ' ( B1 ) 


where a — a/Lo, 


u(/c) = (l + « 2 ) 1 + l/ , (B2) 

and J\(x) is the Bessel function of order 1. For a circu- 
lar disk of area L? the radius is a = Lj\fn. 

Using (Bl) one can easily show that 

f o °°dTC AA (0,T) 


and 




’ JHkoc) ’ 

Jo n 

(z 2 < 1+ w) ) 

= To 

JO K v{k) 


(B3) 


as z — > 0. Prom Eq. (31) it follows that r int (L) vanishes 
for L — *• 0 as 


m tq r(i + v) / M 2M 

intW ~ 4 (1 + 2v)T(\v\)H{v) \2Lo) 


(A4) 


In this limit the product r mt (L)a 2 (L) approaches a con- 
stant : 


lim r int (L)(T 2 (L) ~ 

L — ►O 


7 ot 0 r(i + v) 

2 1 + 2v ‘ 


(A5) 


Appendix B: Effect of Finite Cut-off in 
the Evaluation of the Integral 
Correlation Time 

In this appendix we obtain an expression for the 
fractional error in the integral correlation time fi nt pre- 
dicted by the spectral model for a circular disk of radius 
a as a function of the lag r max at which numerical in- 
tegration of the correlation function is stopped. The 
derivation is a simple extension of the calculations de- 
scribed in Bell and Kundu [2003]. 

The covariance of instantaneous area-averaged rain 
rate averaged over a circular disk of radius a separated 
by a time interval r is 

Caa{ 0,t) = -^^d 2 x^d 2 yc(x-y,r) 

= 2 ’ r ‘ 

x / d 2 ke tk ^ x_y ^c(k,r) . 


If we estimate Ti nt by introducing a finite upper cut-off 
in the r- integral, then the fractional error e(r 0 ;r max ) 
defined as 

e(T 0 ;T max ) = — — — • (B4) 


T int 


is given by the formula 


e ( r 0> r max) 


JT If eX P [-(Tmax/7D)v(«)] 

roodn J?(kq) 

JO K V 2 (k) 


(B5) 

Eq. (B5) is numerically evaluated using Mathematica 
[ Wolfram , 1999] to estimate the correction to the inte- 
gral correlation time discussed in section 3.3. A plot of 
e against the ratio r max /To is shown in Figure Bl for 
typical parameter values. 


Appendix C: Estimate of Random Error 
for the Integral Correlation Time of an 
AR(1) Process 

In this appendix we obtain a simple estimate of the 
random error of the integral correlation time of a nor- 
mally distributed random variable undergoing a station- 
ary, linear, first order autoregressive [AR(1)] process, 
evaluated from a single time series of finite length T. 

Consider the time series {x(t)} of length T generated 
by a single realization of a random variable X(t) under- 
going a linear stationary AR(1) process for which the 
lagged covariance function 

7 (r) = <X'(i + T)X'(i)> 


(Cl) 
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Fractional error as function of cut-off 



Figure Bl. Plot of the fractional error (B5) as function 
of the ratio r max /ro computed for the typical spectral 
model parameters v = —0.25 and Lq = 70 km and three 
values of disk radius a. The curves marked (1), (2) and 
(3) represent the results for the radii a = 1 km, 10 km 
and 100 km. The corresponding true integral correla- 
tion times are 0.052ro, 0.19 tq and 0.65r 0 respectively. 


has the simple exponential form 

7 (r) = a 2 e -|r|/TA . (C2) 


Its estimate, obtained from a finite time series, namely, 



T 

dtx'(t 4- r)x'(t) 


can be thought of as a random variable 

Y{t) = ^ j\t x '{t + r) x '{t) (C3) 

depending on the realization of X(t ). We are interested 
in estimating the random error for the integral correla- 
tion time represented by the random variable 


Tint — 



(C4) 


where the variable Z(r) = Y(r)/Y( 0) represents the 
lagged correlation function. The desired error estimate 
is simply the square root of the variance 


<T 2 [T int ] = (Tj 2 t ) - (Tint) 2 , 


which can then be expressed in the form 

/•Tmajc /*T ma x 

^ Pint] = / / dTMZ'iTjZ'fo)) . (C5) 

Jo Jo 


For simplicity of notation we write 
Y{t) = A,Y(0) = B. 

These random variables can each be expressed as the 
sum of respective mean values and fluctuations: 

A = a -h ck, £ = 6 + /3, 


where 

a = {A)= 7 (r) , b h {B) = cr 2 . 

Up to second order in the fluctuations one finds 

(Z'(ti)Z'(t 2 )) « ^ ^aia 2 + ^-oio 2 ^ , (C6) 

where the subscripts 1 and 2 indicate evaluation at lags 
ri and r 2 respectively. 

Eq. (C6) can be recast in a much simpler form if 
X(t) (and therefore X'(t)) is a gaussian random vari- 
able. Using the relation [see, e.g., Anderson, 1958] 

(x'itox'Mx'wx'Cu)) = 
(x'itOx'ihMx'Mx'fa)) 

+ {X'{t{)X'(t 3 ))(X'{t 2 )X'(U)) 

+ (X'{t x )X'{U)){X'{t 2 )X'{t z )) , 
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which expresses the fourth moment of such a variable in 
terms of products of the second moment, the quantity 

{aia 2 ) = (y'(r 1 )y'(r 2 )) (C7) 

can be reduced to the form 
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